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ABSTRACT 

We investigate seven Monte Carlo algorithms - four old and three new - for construct- 
ing merger histories of dark matter halos using the extended Press-Schechter (EPS) 
formalism based on both the spherical and ellipsoidal collapse models. We compare, 
side-by-side, the algorithms' abilities at reproducing the analytic EPS conditional (or 
progenitor) mass function over a broad range of mass and redshift (z = to 15). Among 
the four old algorithms (Lacey & Cole 1993, Kauffmann & White 1993, Somerville & 
Kolatt 1999, Cole et al 2000), we find that only KW93 produces a progenitor mass 
function that is consistent with the EPS prediction for all look-back redshifts. The 
origins of the discrepancies in the other three algorithms are discussed. Our three 
new algorithms arc designed to generate the correct progenitor mass function at each 
time-step. We show that this is a necessary and sufficient condition for consistency 
with EPS at any look-back time. We illustrate the differences among the three new 
algorithms and KW93 by investigating two other conditional statistics: the mass func- 
tion of the ith most massive progenitors and the mass function for descendants with 
N p progenitors. 
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1 INTRODUCTION 

In the hierarchical structure formation scenario, dark mat- 
ter halos grow by accreting and merging with other halos. 
Statistically modeling halo merger histories is important for 
understanding a diverse spectrum of astrophysical processes 
ranging from galaxy formation, the growth of super-massive 
black holes, to cosmic reionization. 

Numerical simulations aside, the most frequently used 
theoretical framework for studying the build up of dark 



posed for this purpose ( 


see, e.g. , Lacey & Cole|1993| Kauff- 


mann & White 1993; Sheth & Pitman 1997 Sheth & Lem- 


son||1999 Somerville & 


Kolatt|1999| |Cole et al.|2000| |2008| 


Moreno & Sheth||2007 


Neistein & Dekel||2008b). These al- 



matter halos is the Press-Schechter (PS) model (Press & 
Schechter|19 74 ) . This framework is further developed in the 
so-called extended Press-Schechter (EPS) model (Bond et 



gorithms allow one to produce realizations of halo merger 
trees stretching back to high redshifts in a fraction of the 
time that is required for performing and analyzing cosmo- 
logical TV-body simulations of comparable resolution. 

Thus far, most of the commonly used Monte Carlo 
methods are based on the spherical EPS theory. In |Lacey| 
& Cole 1993 (also see |Bond et al. 19911, halo mergers at 



aLlpMlj |Lacey fc Cole|1993"l |Mo fc White|1996||Sheth, Mo 
fc Tormen ||2001| |Sheth fc Tormen|2002| >. For a descendant 
halo of a given mass at redshift zq, the EPS model predicts 
the average mass spectrum of its progenitors at a higher 
redshift z\ (the conditional or progenitor mass function). 

The EPS model provides only statistical information 
about halo mergers and does not specify how progenitor ha- 
los are to be grouped into descendant halos. However, it is 
often useful, particularly in semi-analytic modeling, to have 
actual realizations of the merging history for a large set of 
haloes. A number of Monte Carlo algorithms have been pro- 
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each time step are assumed to be binary: one of the progen- 
itor masses is randomly drawn from the conditional mass 
function, and the other progenitor mass is defined by the dif- 
ference between the descendant halo mass and this first pro- 
genitor mass. Though this seems to be the most natural way 
to generate halo merger histories, it has been pointed out by 
several authors that the binary picture does not reproduce 
the EPS progenitor abundance at earlier times (see, e.g. , 
|Somerville fc Kolatt|1999 1. Moreover, this problem does not 
disappear when the time step is greatly reduced. This fact 
has led to the investigation of alternative Monte Carlo algo- 
rithms with different recipes for building halo merger trees 
in the spherical EPS framework. For example, Somerville & 
Kolatt (1999) find that if the binary assumption is relaxed 
while taking into account the contribution of mass from con- 
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tinuous accretion then the progenitor abundance at large 
look-back times is better reproduced. Cole et al. (2000), on 
the other hand, include diffuse accretion but preserve the 
assumption of binary mergers. More recently, partially due 
to the rapid advances in N-body simulation, various other 
algorithms have been proposed that are either designed to fit 
N-body results (e.g. , |Parkinson, Cole fc Helly ||2008| |Cole 



et al. 2008 Neistein & Dckcl 2008a[) or are based on the 



spherical (Neistein & Dekel 2008b) or ellipsoidal (Moreno 



& Sheth 2007 1 excursion set model. The presence of these 



numerous Monte Carlo algorithms suggests that building a 
Monte Carlo algorithm that is fully consistent with the un- 
derlying EPS model is not unique and can be non-trivial. 
We were motivated to write this paper for a number of 



reasons. First, this is a sequel to our previous work ( Zhang, 



Ma & Fakhouri 2008 1 , which presented an accurate an- 



alytic formula for the conditional mass function for small 
time-steps in the ellipsoidal EPS model. This formula is par- 
ticularly useful as an input for high-resolution Monte Carlo 



simulations of halo merger trees. Earlier formulae (e.g. Sheth 
fcTormen|2"002 | were accurate only for larger look-back red- 



shifts (zi — zo > 0.1). Taking such a large time-step would 
limit the dynamic range in both the progenitor mass and 
redshift that can be covered in a Monte Carlo simulation. 
In addition, until recently, all previous Monte Carlo algo- 
rithms were studied in the framework of the spherical EPS 
model, which is well known to produce inaccurate total (i.e. 
unconditional) halo mass function when compared with sim- 
ulations. This paper will investigate the algorithms in the el- 
lipsoidal model using the formula in Zhang, Ma & Fakhouri 
[1 ( |2008 >. 

Second, as we began to investigate the various Monte 
Carlo algorithms proposed in the literature, we were frus- 
trated by the lack of direct comparison among the different 
methods, each of which has its own range of validity and own 
set of assumptions about how to group progenitors into de- 
scendants (e.g. binary vs multiple progenitors; how the mass 
in progenitors below mass resolution is treated). Moreover, 
it was not always clear why a given algorithm succeeded or 
failed. In this paper, we examine closely the four most fre- 
quently used a lgorithms -|Lacey fc Cole||1993 Kauffrnann| 
|fc White|1993l|Somerville fc Kolatt|1999| [Co!e et al |2000| 



and compare their predictions for the conditional mass func- 
tion over a wide range of progenitor masses and look-back 
redshift (e.g., down to 10 -6 of descendant mass and up to 
redshift 15, much larger than those studied previously). We 
find that only Kauffmann & White ( 1993 1 is fully consistent 



with the EPS model at all look-back time steps. The limita- 
tions and causes of discrepancies in the other three methods 
are discussed. 

Third, in light of the discrepancies in earlier algorithms, 
we investigate three new Monte Carlo algorithms that are all 
constructed to reproduce accurately the EPS predicted con- 
ditional mass function at any look-back redshift. We present 
a consistency criterion that is useful as a general guide for 
building Monte Carlo algorithms: If an algorithm reproduces 
the EPS progenitor mass function for a sequence of sim- 
ulation time-steps between Zi and Zi+i (where i = 0,N), 
then it is guaranteed to reproduce the EPS progenitor mass 
function at any Zj for descendants at any later zu (where 
j, k — 0, iV). This is a necessary and sufficient condition. 

Fourth, the EPS model is an incomplete theory that 



predicts only a subset of statistical properties of halo merg- 
ers. It therefore leaves one with much freedom in how to 
assign progenitors to descendants in a given Monte Carlo 
algorithm. For instance, it is possible to construct differ- 
ent consistent Monte Carlo algorithms that predict different 
statistical merger quantities beyond the conditional mass 
function. Our three new algorithms and KW93 are four ex- 
amples that are degenerate in the conditional mass function 
but are different in other progenitor statistics. In this pa- 
per we illustrate the differences among the models with two 
such statistics: the mass function of the i th most massive 
progenitors and the mass function of progenitors for descen- 
dant halos with N p progenitors. Results from JV-body sim- 
ulations will be needed to constrain these higher-moment 
statistics. Since computing the statistics of progenitor dark 
matter halos in simulations is by itself a major independent 
project, we will focus on the EPS theory and Monte Carlo 
algorithms in this paper and leave the comparison with N- 
body results to a subsequent paper (Zhang, Fakhouri & Ma 
2008, in prep). 

The paper is structured as follows. The EPS formal- 
ism based on both the spherical and ellipsoidal gravitational 
collapse models is reviewed in §!2] In §!3]we discuss three in- 
gredients for how to grow an accurate Monte Carlo merger 
tree: the consistency criterion for reproducing EPS ({ 3.1 1, 



the asymmetry in the EPS progenitor mass function and 
the necessity of non-binary mergers in an algorithm (§3.2), 
and the role of mass resolution and diffuse accretion for pro- 
genitor mass assignment (§3.3). Details of the four old and 
three new algorithms are discussed in §4 and §5, respectively. 
Whenever possible, the resulting progenitor mass functions 
from different algorithms are shown on the same plots for 
ease of comparison. §6 compares the two new progenitor 
statistics that can be used to discriminate among the Monte 
Carlo algorithms that are consistent with EPS. We summa- 
rize our findings in §j7] with a discussion of some recent work 
in this field. 

The calculations in this paper assume a ACDM model 
with fi m = 0.25, D. b = 0.045, h = 0.73, fi A = 0.75, n = 1, 
erg = 0.9. This is the same cosmology used in the Millennium 
simulation ( Springel et al.|2005 l. 



2 AN OVERVIEW OF EPS 

In this section we present a brief overview of the EPS the- 
ory based on both the spherical and ellipsoidal gravitational 
collapse models. We often refer to the two models in par- 
allel as the spherical and ellipsoidal EPS models, with the 
understanding that the ellipsoidal version is based on the 



excursion set formalism of Bond et al. ( 1991 1. The emphasis 



here is on the conditional mass function, which is the main 
statistical quantity used to generate progenitors in merger 
tree algorithms. For a more complete and pedagogical review 



of EPS, see Zentner (20071 and references therein 



2.1 EPS Based on the Spherical Collapse Model 

The Press-Schechter (PS) model provides a framework for 
identifying virialized dark matter halos. It is assumed that 
the seed density perturbations that grow to form these ha- 
los are characterized by an initially Gaussian random den- 
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sity field with larger fluctuations on smaller spatial scales. 
This latter assumption allows one to use S(R) = a 2 (R), the 
variance of the linear density fluctuations^] smoothed over 
spatial scale R, as a proxy for the spatial scale R. Moreover, 
since a given spatial scale is related to a unique mass scale 
M(R) via the mean density of the universe p, one can use 
R, M , and S interchangeably as measures of scale. 

The density field smoothed over a given scale M is 
given by 5 m — pu/p — 1 where pu is the average den- 
sity within the smoothing scale R. In the EPS model, the 
linear density field centered at a given point in the initial 
Lagrangian space traces out a random walk (referring to 
a Markovian process^] as the smoothing scale is reduced. 
Starting from a large smoothing scale, a virialized dark mat- 
ter halo is assumed to form at the given spatial coordinate 
when the linear 8 m crosses a critical value for the first time; 
the mass of the halo is determined by the smoothing scale 
at first-crossing. In the spherical EPS model, the critical 
over-density is given by the spherical collapse model and is 
a constant S c = 1.69 independent of mass scale. 

In the above description, as a result of gravitational in- 
stability, the density field grows with time as a linear func- 
tion of its initial value, i.e. , 5m(z) — 8m(0)D(z), where 
D(z) is the standard cosmology-dependent linear growth 
factor satisfying D(z = 0) = 1. In practice, one usually fixes 
the value of 5m at some reference time (e.g. today: 8m(0)) 
and evolves the critical over-density to identify virialized 
halos at earlier redshifts. We denote this time-dependent 
critical over-density by u(z) = S c /D(z). Note that a lower 
redshift corresponds to a smaller u)(z), implying that larger 
halos form at later times, in accordance with the hierarchical 
structure formation scenario. 

Under the assumption of Gaussian statistics, the EPS 
framework allows one to compute the first crossing distri- 
bution f(S(Mi), zi\S(Mq), zo). Of the set of random walks 
that begin at 5m = w(zo), the first crossing distribution is 
the fraction of these random walks that first cross the crit- 
ical over-density uj(zi) at scale S(Mi), where Zi > zo and 
g(Mi) > S{M ) (i.e. Mi < M ). It can be shown (|Lacey fc 



|Cole|1993| > that the first crossing distribution in the spheri- 
cal EPS model has the form 



f(S(M 1 ),z 1 \S(M ),z )dAS 



(1) 



2^ AS 3 / 2 



exp 



(Au 



2AS 



dAS 



where Auj = u(zi) - lu(z ) and AS = S(Mi) - S(M ). 

The first crossing distribution can be reinterpreted as 
the conditional mass function P(Mi, z\\Mq, zq), which is de- 
fined to be the mass fraction of a descendant halo of mass 
Mo at redshift zo that originates from a progenitor halo of 
mass Mi at redshift z\: 

P{M 1 ,zi\M ,z )dM 1 = -f(S(Mi),zi\S(M ),zo)dAS (2) 

Note, in particular, that P(Mi, z±\Mo, zq) is the mass- 
weighted conditional mass function as it represents the merg- 



1 In this paper, the variance of the density fluctuation is calcu- 
lated using the fitting formula of the linear mass power spectrum 
from |Eisenstein fc Hu|1998| 

2 Strictly speaking, this is only true when the smoothing window 
function is a top-hat in Fourier space. 



ing history of a unit of mass. The average number of progen- 
itors of mass Mi at zi associated with the formation of each 
descendant halo of mass Mo at zo is given by the number- 
weighted conditional mass function (j>(Mi, zi\Mo, zo), which 
is simply related to the mass-weighted conditional mass 
function P(Mi, zi\M , z ) by 

M 



<t>(Mi,zi\Mo,zo) 



Mi 



P(Mi,zi\M ,zo) 



(3) 



For brevity, we often refer to the number-weighted condi- 
tional mass function <f)(Mi, zi\Mo, zo) as the progenitor mass 
function, and denote it simply as 4>(Mi\Mo) with z and zi 
specified elsewhere in paper. This quantity is sometimes de- 
noted as dN(Mi, zi\Mo, zo) / dMi in the literature. 

2.2 EPS Based on the Ellipsoidal Collapse Model 

The original Press-Schechter theory was based on the spher- 
ical collapse model. The unconditional mass function in this 
model is well known to have an excess of small halos and 
a deficit of massive halos in comparison with simulation re 
suits (e.g. ,|Lacey fc Cole||i"994l |Gelb fc Bertschinger|[l99^ 
Ma fc Bertschinger||1994[ |Tormen||1998| |Sheth fc Tormenl 



on 



19991). This discrepancy arises because halo collapses are 
generally triaxial rather than spherical ( Doroshkevich|1970 



Bardeen et aLl|l986l [Sheth, Mo fc Tormen ||2001[ |Sheth fc| 



Tormen||2002 1. In the spherical collapse picture, the virial 



ization of a dark matter halo is purely determined by the 
density-contrast on the scale of the halo mass. This assump- 
tion is too simplistic because dark matter halos generally 
have non-zero ellipticity and prolateness, and the condition 
for virialization should be determined by both the density- 
contrast and the halo shape parameters. 

By assuming that a dark matter halo virializes when its 
third axis collapses, ISheth, Mo & Tormen ( 20011) find a new 



criterion for virialization that depends on the ellipticity and 
prolateness of a dark matter halo in addition to its density 
contrast. In practice, this condition can be simplified either 
by averaging over its dependence on the shape parameters, 
or by fixing the shape parameters at their most likely values 
for a given over-density. By doing so, these authors obtain a 
fitting formula for the scale- dependent critical over-density, 
or barrier, in contrast to the scale-independent 8 C of the 
spherical collapse model. It is parameterized as (Sheth & 
Tormen|200"2"l ): 

5f\S(M), Z \ = ^ 1 5 c \l + 0( ri u)- c '} (4) 

where q = 0.615, = 0.485, 7 = 0.75, v = cj 2 (z) / S '(M) , and 
M is the halo mass. In this ellipsoidal collapse model, the 
scale-dependence is such that the formation of small halos 
is delayed, thereby reducing their abundance and providing 
closer agreement with the unconditional mass function in 
simulations than the spherical model. 

To compute the conditional mass function in the ellip- 
soidal EPS model, one would need the equivalent of the first- 
crossing distribution eq. |TJ. The exact analytical form of 
eq. |TJ , unfortunately, is valid only for the scale-independent 
constant barrier 8 C of the spherical EPS model. |Sheth fc] 



Tormen (2002) have presented a Taylor-series-like approxi- 



mation for the ellipsoidal model, but Zhang, Ma & Fakhouri 
[1 1 |2008[ ) show that this form works well for large zi — zo but is 
invalid for small zi — zo- As the construction of an accurate 
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ellipsoidal Monte Carlo merger tree algorithm requires ac- 
curate knowledge of the ellipsoidal progenitor mass function 
at small time-steps, it is crucial that this matter be resolved. 
This was done inlZhang, Ma & Fakhouri I (120081). Using 



the scale-dependent critical over-density of Sheth & Tormen 



| 2002[) and the techniqu e of |Zhang fe Hui| (|2006f , |Zhang, 
|Ma fe Fakhouri~| ( |2008[ ) derived an accurate form for the 
progenitor mass function of ellipsoidal EPS model for small 
time steps (Az < 0.1), which can be written as: 



0(Mi|M o ) = 



M dS(Mi) AqAlu 



Mi dM 1 ASV^AS 

2^S ( A ° AlJ + Al ^^)' 



(5) 



cxp 



+ A 2 S 3/2 exp 



where A = 0.866(1 - 0.133i>o 



Ai 



r(3/2) 



0.308^, 



o 



0.0373^; 



f'o 



ui 2 (z )/S(M ), and g 



AS/S(M ) . Not e that unlike eq. (15) of |Zhang, Ma fe 
Fakhouri (20081, we have not neglected the small AqAui 



term in the exponent because it is important for trac- 
ing the massive progenitors (small AS). Two other fea- 
tures are worth noting. First, unlike in the spherical EPS 
model, <j)(M\\Mo) in eq. |5| depends weakly on the red- 
shift 20. Second, due to the intersections of barriers at the 
low mass enoj^] eq. ([HJl turns unphysical (i.e. , Ao < 0) 
when S(Mo) ^ 30uj 2 (zq), i.e. , when the descendant mass 
is much smaller than the typical halo mass at zo. In our 
Monte Carlo simulations discussed below, whenever the sec- 
ond feature becomes a problem (which occurs very rarely), 
we do not generate any progenitors for the halo in the next 
time step. As we will show in Sj5] this procedure only mildly 
affects the progenitor abundance at the very low mass end. 
Eq. (|5| provides a closer match to the merger rates deter- 
mined from iV-body simulations ( Zhang, Ma fc Fakhouri 
2008[ ), but the agreement was not perfect, perhaps due to 
the non-Markovian nature of numerical simulations. 



3 INGREDIENTS FOR GROWING HEALTHY 
MONTE CARLO MERGER TREES 

As discussed in the introduction, the EPS model only pro- 
vides a subset of statistical information about dark matter 
halo merger histories. For example, the EPS progenitor mass 
function <p(M\Mo) (eq. [3] for spherical and eq. [5] for ellip- 
soidal) gives the average mass spectrum of the progenitors 
for the descendant halos. However, it is often useful, espe- 
cially in semi-analytical modeling, to have an actual Monte 
Carlo realization of the formation history for a large set of 
halos. Of particular interest is the merger tree of individ- 
ual halos, which provides the hierarchical links among the 
progenitors and their descendants. Since the EPS model it- 
self does not specify explicitly how to group progenitors into 
descendants, in each time-step in a Monte Carlo algorithm, 
assumptions must be made about the number of progeni- 
tors and their mass distributions to be assigned to a given 
descendant. 



The earlier Monte Carlo algorithms (e.g., 


Lacey & Cole 


1993| IKauffmann & White|1993| |Somerville & Kolatt|1999 


Cole et al.|2000 


i for merger tree constructions share a simi- 



lar overall structure: A descendant halo of mass Mo at some 
redshift zo (typically zq = 0) is chosen. The EPS progen- 
itor mass function, <j>(M\Mo), is then used to generate a 
set of progenitors at some earlier redshift, using the rules of 
the given algorithm. In the next time-step, these progenitors 
become descendants, and each is assigned its own set of pro- 
genitors at an earlier redshift using cj)(M\Mo). This process 
is repeated out to some early redshift and for a (typically 
large) number of halos of mass Mo at the starting zo- 

The existence of a number of diverse Monte Carlo algo- 
rithms (see further discussion in Q in the literature implies 
that the above process is, in fact, not unique and can be 
quite subtle. We now explore some of these subtleties and 
the key ingredients for constructing a healthy merger tree. 



3.1 A Criterion for Consistently Reproducing the 
EPS Progenitor Mass Function 

We consider a Monte Carlo algorithm to be consistent with 
EPS if the merger trees it produces can reproduce the EPS 
progenitor mass function (j>(Mi, Zi \Mq, zo) exactly for any 
set of {Mi, zi, Mo, zo} regardless of the number or size of 
the simulation time-steps between zo and Z\. 

Clearly, to be consistent with EPS, a Monte Carlo 
algorithm must necessarily reproduce the EPS-predicted 
<f>(M\Mo) exactly at adjacent time steps. We now show that 
this is also a sufficient condition for the Monte Carlo method 
to reproduce cj>(M\Mo) exactly at any look-back time re- 
gardless of the number, or width, of intervening time-steps. 
This condition is important because it simplifies the analysis 
of Monte Carlo algorithms: the failure of a given algorithm to 
reproduce faithfully the EPS <j>(M\Mo) at a particular red- 
shift or mass range necessarily implies that the algorithm 
fails to reproduce the progenitor mass function (in either 
amplitude or shape or both) across a single time step. 

We start with the first crossing distribution eq. ([!]) and 
note that due to the continuous nature of the random walk, 
it obeys the following identity at different look-back times: 

f(S(M),z\S(M ),z ) (6) 

S(M) 

dS'f(S(M), z\S(M'), z') f(S(M'),z'\S(M ), z ) 

S(M ) 

for any zo < z < z. This relationship is true in both spher- 
ical and ellipsoidal EPS models because both variants are 
based on barrier crossings of random walks. Note that in the 
ellipsoidal model, eq. |6| is a property of only the exact first- 
crossing distribution, which is well represented by eq. |5| for 
small look-back times but not the Taylor-series-like approx- 
imation of Sheth & Tormen (20021. Also note that eq. |6| 



may not be strictly satisfied at the very low mass end due 
to the intersections of barriers in the ellipsoidal model. As 
we will show in [J5] this only causes a minor problem on very 
small mass scales. 

Using eqs. |2| and |3| to relate / to the progenitor mass 
function <t>, we then obtain 



see appendix A of Sheth fc Tormen|2002 for more details 



(j>(M,z\Mo,z ) 



(7) 
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= / dM'<j)(M,z\M',z')(f>(M',z'\M ,z ) . 

J M 

Setting z = zq + Az and z = zo + 2Az, we see that eq. |7| 
implies that if a Monte Carlo method generates progeni- 
tors exactly according to the progenitor mass function of 
EPS at each time step Az, then the Monte Carlo progeni- 
tor mass function should agree with the EPS prediction at 
any look-back time 2 — 2:0. We stress that cj>[M [Mo) must 
be reproduced exactly, that is, in both the overall shape 
and normalization of <p(M\Mo). This consistency condition 
is both necessary and sufficient. 

An additional feature to note is that consistency is pos- 
sible in the presence of a mass resolution limit M res (dis- 
cussed further in j |3.3[ ). Eq. |7| shows that 4>(M, z\M ,z ) 
does not depend on masses outside of the range [M, Mo]. 
Thus if a Monte Carlo algorithm reproduces 0(M|Mo) for 
all M > M rcs in single time-steps, it will consistently repro- 
duce </>(M|M ) at M > M res for any 2 — 20. 



3.2 The Asymmetry of EPS and Binary Mergers 

The simplest way to group progenitors into descendants in 
a Monte Carlo algorithm is through binary mergers, i.e. , 
each descendant halo of mass Mo is composed of two pro- 
genitors of mass Mi and Mo — Mi . This assumption is used 
in, e.g., Lacey & Cole ( 19931. This simple scenario, however, 



will necessarily fail to reproduce both the spherical and el- 
lipsoidal EPS progenitor mass functions. This is because if 
all descendants were the products of binary mergers, then 
0(Af|Mo) would be symmetric about Mo/2 for infinitesimal 
Az. This is simply not the case in EPS. 

We illustrate the asymmetry of the EPS <^(M|M ) in 
Fig. [I] for a descendant halo of mass 1O 13 M0 at 20 = and 
a look-back time of 2 = 0.02 (which is the typical time-step 
used in our Monte Carlo simulations; see Sec. 4). The solid 
black curve shows the total (f>(M\Mo), while the red dashed 
curve shows the symmetric part S ym (M|M ) defined by 



4>{M\M Q ) = <^ ym (M|M ) + ^ aS ym(M|M ) 



(8) 



where the left side (M < M /2) of ^ sym (M|M ) is defined 
to be identical to <fi(M\Mo) and the right side is defined to 
be simply the reflection of the left half about the mid point 
Mo/2. The second term </> aS ym(M|Mo) is then the residual 
of (f> after subtracting out <f> sym . The figure illustrates that 
it is not possible for all progenitors with M > Mo/2 to 
have binary-paired progenitors of mass Mo — M < Mo/2. In 
particular, we find that for sufficiently small look-back times 
(e.g. 2 = 0.02 used in Fig. [TJ, 0(M|M o ) > <^ sym (M|M ) 
when M /2 sC M <, 0.97M and <f>(M\M ) < <^ sym (M|M ) 
when M > 0.97 M (see the pop-up in Fig. [I}. That is, there 
are slightly fewer progenitors with masses below Mo/2 than 
above, except near the end points (below 0.03Mo and above 
0.97Mo) where the trend is flipped. 

Even though the asymmetry is typically small ( <^ aS y m < 
0.1</> S ym out to Mo — M res ), an accurate algorithm must in- 
clude non-binary progenitor events. These can be descen- 
dants with either a single progenitor or multiple (N p > 2) 
progenitors, as will be seen in the new algorithms discussed 
in f|5] below. This fact was emphasized by |Neistein fc Dekel| 
(2008bJ. These authors construct a mass conserving consis- 
tent Monte Carlo algorithm that produces a large number 




0.4 0.5 0.6 

M/M„ (A/„ = 10 13 Af o ) 

Figure 1. An illustration of the asymmetry in the 
number-weighted conditional (or progenitor) mass function 
<j>(M, z\Mo, zo) of the spherical EPS model for a descendant halo 
of mass Mq = 1O 13 M0 at zo = and a look-back redshift of 
z — zo = 0.02. The red dashed curve shows the symmetric part of 
4>(M\Mo), <fisym(M\Mo), whose right side is simply the reflection 
of the left side. The figure indicates that some progenitors of 
masses larger than Mq /2 do not have companions in the simplest 
binary scheme. The pop-up is a zoom-in on the right-most part 
of the plot and illustrates that the red dashed curve exceeds the 
solid curve at M 0.977M . 



of non-binary descendants. However, one intuitively expects 
that more mergers will be binary as 21 — 20 — > 0. This in- 
tuition is supported by results from the Millennium simu- 
lation ( |Fakhouri fe Ma 2008), which show that the binary 
assumption becomes increasingly valid down to smaller M res 
as z\ — zo is made smaller. This result suggests that the 
Markovian nature of the standard EPS model with a top- 
hat smoothing window may need to be modified to account 
for the correlated sequences of mergers occurring in simula- 
tions ( jNeistein fc Dekel|2008b"l |Zentner|2007l ). 

3.3 Mass Resolution, Diffuse Accretion, and Mass 
Conservation in Monte Carlo Algorithms 

In the EPS model, all the mass in the universe is assumed 
to be in dark matter halo^] Although the mass-integral of 
the (unconditional) mass function in this model is finite, 
the number-integral is unbounded; that is, EPS predicts a 
preponderance of very low mass halos. Thus, any practi- 
cal Monte Carlo algorithm must necessarily assume a lower 
mass cutoff, the mass resolution M rcs . 

For a nonzero M rea , a halo's merger history at each time 
step can be thought of consisting of mass in the form of 
resolvable progenitor halos and a reservoir of mass due to 
"diffuse" accretion that is the aggregate contribution from 



4 This is not exactly true in the ellipsoidal EPS model. See ap- 
pendix A of |Sheth fc Tormen|20"02| 
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all sub-resolution progenitors. This technical distinction is 
introduced for ease of implementing the Monte Carlo meth- 
ods. It will, however, play a more physical role when we 
compare the results with TV-body simulations, which has its 
own mass resolution as well as a possibly physical diffuse 
component consisting of tidally stripped dark matter parti- 
cles. In this paper we use AM to denote this diffuse accretion 
component, which we define to be 



AM = Mo - Mi 



(9) 



where Mi are the masses of the progenitors above M Tes and 
Mo is the mass of the descendant. 

We call a Monte Carlo algorithm mass conserving if 
each descendant and its progenitors produced by the algo- 
rithm satisfies ^ - Mi ^ Mo . Monte Carlo algorithms are 
generally expected to be mass conserving, but we note that 
this is not a necessary condition for reproducing the EPS 
progenitor mass function because the latter is a statistical 
measure of merger properties. In two of our new algorithms 
below (methods A and B in Sj5|, a small fraction of the de- 
scendants can have Mi > Mo- We allow this to simplify 
the description and implementation of our algorithms. We 
have experimented with redistributing these excess progen- 
itors among other descendant halos in a mass-conserving 
manner and found it not to modify significantly the re- 
sulting merger statistics. In addition, it may appear that 
y\ Mi > Mo is unphysical. We have found, however, that 
a non-negligible fraction of halos in iV-body simulations in 
fact have AM < 0, perhaps as a result of tidal stripping. 
This point will be discussed in greater detail in our next 
paper. 

We note that for a Monte Carlo algorithm that is consis- 
tent with EPS, the mean value of AM per descendant halo 
of mass Mo (i.e., averaged over all descendants in a given 
time-step) is, by construction, related to the mass resolution 
by 



(AM) 



M(f>(M\Ma)dM . 



(10) 



For a given 4>(M\Mo) and M res , (AM) is therefore specified. 
The distribution of AM, however, can differ greatly among 
different algorithms; that is, there is much freedom in how to 
assign the amount of diffuse accretion to individual descen- 



dants in a given time-step. For instance, Cole et al. (20001 



assumes a delta-function distribution for AM (see j ]4.4| for 
details), while most of other methods, including our new 
methods discussed in fj5] have broader distributions. 



4 COMPARISON OF FOUR PREVIOUS 
MONTE CARLO ALGORITHMS 

In this section we examine four existing Monte Carlo al- 
gorithms for generating merger trees: |Lacey fc Cole| (|1993[) 



Kauffmann fc White| (p93| [ KW93 ] , |Somerville fc 
1999| l [SK99], and |Cole eFaLlpOOOl ) [C00]. This set 



[LC93] , 
|Kolatt 

is by no means complete, but these are four of the most fre- 
quently used algorithms in the literature. The purpose here 
is to compare these well known algorithms side-by-side and 
to illustrate the mass and redshift ranges for which each 
method succeeds and fails in matching the spherical EPS 



model. This not only benefits the current users of the meth- 
ods, but also prepares us for incorporating the ellipsoidal 
EPS model into the successful method (KW93), which will 
be compared with our new methods in Sec. 5. 

We review each algorithm in a subsection below and 
compare the resulting progenitor mass functions <^>(M|Mo) 
with the spherical EPS prediction for look-back redshifts 
ranging from 0.24 to 15. In Figs. |2|4~] we plot the progenitor 
mass functions produced by all four methods, along with the 
analytical EPS prediction, on log-log plots for three descen- 
dant masses (10 12 , 10 13 , 1Q 14 Mq) and four look-back times 
(zi — zo = 0.24,2.07,7 and 15). To ease comparison, we 
also plot the ratio between each Monte Carlo result and 
the EPS prediction on a linear- log plot. As Figs. |2|4"| clearly 
show, of the four algorithms, only KW93 is able to match 
the spherical EPS 4>(M\Mo) for all z — zq. We will explore 
why each algorithm fails below and discuss the care that 
must be taken when implementing KW93. A summary of 
the four algorithms, their discrepancies, and the causes of 
the discrepancies is given in Table [l] 

In our Monte Carlo simulations, we generally keep track 
of all progenitors down to O.OOlMo at each time step for 
a descendant halo of mass Mo. This large dynamic range 
allows us to predict reliably the progenitor abundance even 
for a very large look-back time (zi — zq ~ 15). To speed 
up the algorithm, we take each time step to be a constant 
difference in the barrier height Auj(z) = uj(z + Az) — u(z) 
(where Auj ~ Az at low z), which is chosen to be about 0.02 
for LC93, KW93, SK99, and 0.003 for COO at z = 0. The 
progenitor mass function of a given descendant halo mass is 
then identical for each time step and does not have to be 
recomputed. Numerical convergence is tested by changing 
the time-steps used in the simulation: our results do not 
change. 



4.1 Lacey & Cole (1993) 

The algorithm proposed by LC93 makes two important as- 
sumptions: all mergers are binary (before mass resolution 
is imposed), and the descendant mass Mo is the sum of 
the two progenitor masses Mi and M2 (where Mi ^ M2 
in our convention). For each small look-back time step and 
for each descendant, a progenitor mass is randomly chosen 
according to the mass-weighted conditional mass function 
cq. Q, and the mass of the other progenitor (which can be 
larger or smaller) is simply set to be the difference between 
Mo and the first chosen progenitor mass. If the less massive 
progenitor M2 falls below a chosen mass resolution M ICS , or 
equivalently, Mi > M - M Ica , then Mi is kept but M 2 , 
being a sub-resolution progenitor, is discarded. This results 
in single-progenitor halos which we label as "1 — > 1" events. 
In this notation, binary mergers are "2 — > 1" events. When 
a smaller time-step is used in LC93, the ratio of 2 — > 1 to 
1 — ► 1 events decreases. 

We find that random progenitor masses can be easily 
generated using the parameter transformation: 



erf 



[Alu/ v / 2[S(M 1 )-S(M )]} 



(11) 



The parameter x has a uniform probability distribution be- 
tween and 1 and can be quickly generated using any 
random-number generator. A simple inversion then yields 
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Algorithm 


Overview 


Discrepancy in progenitor mass function <p(M \Mq) 


Reasons for Discrepancies 


LC93 


Binary and 1-to-l 
AM < M rcs 


Overestimates </>(M\Mq) by large factors when the 
look-back time is large, i.e. , z\ — zq > 1 


Binary assumption fails to reproduce EPS 
4>(M\Mq) asymmetry. 


KW93 


Multiple mergers 
AM ^ 


None 




SK99 


Multiple mergers 
AM ^ (can be 
bigger or smaller 
than M res ) 


Typically over-predicts the abundances of small 
progenitors ( < 10% of the descendant halo mass) 
by a factor of ~ 2 for zi — zq 1. This discrepancy 
propagates to smaller mass scales for larger 
look-back times. 


Truncation of 4>{M\Mq) fails to reproduce 
its shape exactly. 


COO 


Binary and 1-to-l 
AM is a constant 
given by equation 

& 


Works reasonably well for a large range of the 
look-back time but significantly underestimates 
<f>(M\Mo) at high mass ends, particularly when the 
look back time is large (zi — zg S> !)■ 


Binary assumption fails to capture 
asymmetry of EPS </>(M|M ); fixed AM 
yields 1-to-l events that do not accurately 
reproduce the high mass end of (j>(M\Mo). 



Table 1. A scorecard for the four old Monte Carlo algorithms discussed in [jl] We note that the 1-to-l events in LC93 and COO are 
actually binary mergers involving a secondary progenitor with mass below Af rcs . Since these progenitors are below the resolution limit 
they are not counted as progenitors but as diffuse mass AM. 




M/Mo (M = 1O 12 M ) 

Figure 2. Comparison of the progenitor (or conditional) mass functions <f>(M, z\Mo, Zo) that we generated using the four previous Monte 
Carlo algorithms by LC93 (red solid), KW93 (orange dot-dashed), SK99 (blue dashed), and COO (green dotted), and the predictions 
of the analytic spherical EPS model (black solid). The four panels show four look-back redshifts [z — zq = 0.24,2.07,7 and 15) for a 
descendant halo of Mo = 1O 12 M0 at zo = 0. For clarity, we plot in the sub-panel below each panel the ratios of the Monte Carlo result 
and the EPS prediction. One can see that KW93 is the only accurate algorithm for all z. Note that different ranges of M/Mo are shown 
in each panel since the progenitors have progressively smaller masses at higher redshifts. 



progenitors distributed according to the mass- weighted con- 
ditional mass function. 

The red solid histograms and curves in Fig. [2] -[^com- 
pare the progenitor mass functions generated using the LC93 
algorithm with the predictions of the spherical EPS model 
(solid black curves). For all three descendant halo masses 
shown (10 12 , 10 13 and 1O 14 M0), we see close agreement for 
small look-back times such as z\ — zq — 0.24, but LC93 pro- 
duces an excess of progenitors at larger look-back times, and 
the discrepancy worsens, reaching an order of magnitude by 
z\ — Zo = 15. We believe this discrepancy is due to the binary 
nature of LC93: the number of progenitors with mass M is 



equal to the number of binary companions of mass Mq — M. 
Thus the LC93 Monte Carlo algorithm generates a progen- 
itor mass function after one time step that is symmetric in 
the left and right sides, which will not match the asymmet- 



ric nature of the EPS cf>(M\Mo) discussed in Sec. 3.2 and 
shown in Fig. [T] This discrepancy is amplified after many 
time-steps when the look-back time becomes large. 

Finally, we note that the authors of LC93 also consider 
another way of drawing the first progenitor mass from the 
mass-weighted conditional mass function, which is to draw 
it from the mass range of [Mo/2, Mo] instead of [0,Mo]. In 
practice, we find that this slightly modified version of LC93 
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Figure 3. Same as Fig. [2] but for a descendant halo of 10 13 Mq. 
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M/M (M = 10 14 M o ) 

Figure 4. Same as Fig. [2] but for a descendant halo of 10 14 M(; 



generates very similar results, and our above discussion is 
valid. 



4.2 Kauffmann & White (1993) 

For each time-step in the KW93 algorithm, a large number of 
progenitors are generated across many progenitor mass bins 
for a fixed number of descendant halos of the same mass. The 
number of progenitors in each mass bin is determined by the 
progenitor mass function of the descendant halo mass, and 
rounded to the nearest integer value. These progenitors are 
then assigned to the descendant halos in order of decreasing 
progenitor mass. The target descendant halo is chosen with a 



probability proportional to its available mass (i.e. the mass 
not yet occupied by progenitors), and with the restriction 
that the total mass of the progenitors in a descendant halo 
cannot exceed the descendant mass. This procedure allows 
one to work out all the merger configurations and their fre- 
quencies for one time step and for different descendant halo 
masses. This information is then stored and used repeatedly 
for determining the progenitors of a halo at each time step. 

To speed up the implementation of KW93, we divide the 
look-back time into steps with equal spacing in the barrier 
height Auj as discussed earlier. The progenitor mass function 
for a fixed descendant halo mass is then identical for every 
time step and only has to be calculated once. We store the 
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ensemble of progenitors and their merger configurations for 
each descendant halo mass bin. In a Monte Carlo simulation, 
we randomly select one merger configuration from the many 
stored ones for a descendant halo at each time step. 

In practice, we find that extreme care must be taken 
to avoid numerical problems in KW93. First of all, this al- 
gorithm requires a large number of progenitor mass bins in 
the neighborhood of Mo because <f>(M\Mo) is sharply peaked 
near Mi ~ Mo for small time-steps. Interestingly, we find 
that if the mass range of [M res ,Mo] is simply divided into 
evenly-spaced logarithmic bins, this method is not accurate 
even when the number of mass bins is as large as 2000, 
which already requires more than ~ 50000 descendant halos 
to guarantee that the integer rounding does not introduce 
a significant error to the progenitor number in each bin. As 
a result, a large amount of computer memory is necessary 
to repeat this procedure for descendant halos of different 
masses. The improved mass bin configuration that we end 
up using will be introduced in S|5] Using that setup, we find 
that only 200 bins are required to reproduce accurately the 
EPS progenitor mass function over large look-back times. 

The second problem is that KW93's scheme for assign- 
ing progenitors to descendant halos is somewhat ambiguous 
and does not guarantee that all the progenitors can be as- 
signed. Fortunately, we find that this problem usually does 
not arise when the ensemble of progenitors is large. For each 
descendant halo mass, we use ~ 8000 descendant halos to 
determine the merger configurations of the progenitors. 

The orange dash-dotted curves in Fig. [2] - [4] compare 
the progenitor mass functions generated using the KW93 
algorithm with the predictions of the spherical EPS model 
(black). The results show very good agreement. Since KW93 
reproduces the exact EPS progenitor mass function at every 
time-step, it is expected to be consistent with EPS at any 
Zi — 20 according to the discussion in §3.1| 

4.3 Somerville & Kolatt (1999) 

Somerville & Kolatt (1999) [SK99] point out that the as- 
sumptions of binary mergers and Mo = Mi + M2 made in 
LC93 lead to an overestimate of the progenitor abundance at 
high redshift. They first attempt to remedy this problem by 
preserving the binary assumption while allowing the mass 
below the resolution limit M ros to be counted as diffusely 
accreted mass AM (see £ 3.3 ) . They show, however, that 



this "binary tree with accretion" method fails in the oppo- 
site direction, under-producing the progenitor mass function 
relative to the spherical EPS prediction. This discrepancy 
arises partly because whenever two progenitors are chosen 
in this method, the remaining mass is assigned to AM re- 
gardless of whether it is above or below M res . Thus the EPS 
<f)(M\Mo) is not faithfully reproduced: the binary tree with 
accretion method yields an excess of accreted mass and a 
corresponding shortage of low-mass halos. 

SK99 then consider a natural extension of this method, 
in which both assumptions made in LC93 are relaxed. In 
this "N-branch tree with accretion" algorithm, each descen- 
dant halo is allowed to have more than two progenitors for 
every simulation time-step. To guarantee that the total mass 
of the progenitors does not exceed that of the descendant, 
each subsequent progenitor mass is randomly chosen from 
the mass-weighted conditional mass function truncated to 



the maximally possible progenitor mass. This procedure is 
repeated until the descendant halo cannot contain any more 
progenitors with masses above M rea , and the remaining mass 
deficit is assigned to diffuse accretion AM. 



The parameter transformation of eq. (Ill is also ap- 



plicable for SK99. The probability distribution of x is still 
uniform, but the upper limit of x can now take on any value 
between and 1 depending on where the conditional mass 
function is truncated. 

The blue dashed curves in Fig. [2] - [4] compare the pro- 
genitor mass functions generated using the N-branch tree 
algorithm of SK99 with the predictions of the spherical EPS 
model (black). It is interesting to note that the sign of the 
discrepancy is now opposite to that of LC93: SK99 produces 
an excess of low-mass ( <, O.lMo) progenitors by up to a fac- 
tor of ~ 2 for small look-back times, but it does a better job 
than LC93 at high redshifts. However, it is noteworthy that 
even at high redshifts, discrepancies of up to a factor of ~ 2 
are still present for small progenitor masses. 

We believe that the use of a truncated progenitor mass 
function in SK99 is at least a partial cause for the over- 
prediction of small progenitors. Since the distribution of 
progenitors (in particular, the upper limit for the progenitor 
mass) depends on the sum of the masses of the progenitors 
already picked out for the current halo, the order in which 
progenitor halos are randomly pulled out matters in this 
method. Halos more massive than the truncation limit are 
effectively discarded instead of being randomly selected and 
placed in, for example, new descendant halos. This proce- 
dure tends to preferentially skew the progenitor mass func- 
tion at small time steps towards more low mass progenitors 
and fewer high mass progenitors. 



4.4 Cole et al. (2000) 



Similar to SK99, |Cole et al.| pOOO] ) [COO] treats the mass 
in progenitors smaller than the mass resolution M res in the 
Monte Carlo simulation as accreted mass, but unlike the 
N-branch tree model in SK99, only a maximum of two pro- 
genitors are allowed per descendant. The amount of accreted 
mass gained in one time-step, AM, is fixed to a single value 
and is calculated by integrating the mass-weighted condi- 
tional mass function from to M rcs : 



AM = 



A -fro 



M<f>(M\M Q )dM , 



(12) 



where Mo is the descendant mass. The progenitors are drawn 
from the lower half of the progenitor mass function between 
M res and Mo/2 according to the average number of progen- 
itors in that range: 



M /2 



<t>{M\Ma)dM. 



(13) 



The simulation time-step is chosen to be small enough so 
that p <C I (note that it is for this reason that we use Az — 
0.003 when implementing COO). 

The COO merger tree is generated with the following 
steps: A random number x between and I determines 
whether a descendant halo has one progenitor (if x > p) or 
two progenitors (if x ^ p). In the case of a single progenitor, 
its mass is Mi = Mo — AM. In the case of two progenitors, 



© 2006 RAS, MNRAS 000, [T]-[20l 



10 Jun Zhang, Onsi Fakhouri, and Chung-Pei Ma 



the mass of the smaller progenitor, M2, is chosen randomly 
between M res and Mo/2 according to the progenitor mass 
function. The larger progenitor is then assigned a mass of 
Mi = Mo — M2 — AM. Since p <C 1, most descendants form 
via 1 — > 1 events rather than 2 — » 1 events. To improve the 
speed of this algorithm, we precompute and store the binary 
merger rates and diffuse accretion mass fractions for a single 
time step for different descendant mass bins. 

The green dotted curves in Fig. [2]- [4] compare the pro- 
genitor mass functions generated using the COO algorithm 
with the predictions of the spherical EPS model (black). 
The agreement is noticeably better than LC93 and SK99. 
The largest discrepancy occurs at the high mass end at large 
z\ — zo, where COO under-predicts the progenitor number at 
z\ by more than a factor of two for group-to-cluster size 
descendants at z with Mo > 1O 13 M . 

At least two problems contribute to this discrepancy: 

(i) Since AM is fixed to one value (eq. 12 1, the mass of 
the progenitor for 1 — > 1 descendants is also a fixed value: 
Mi = Mo — AM. This is an over-simplification that com- 
presses the high mass end of 0(M|Mo) into a delta function. 

(ii) For descendants with binary progenitors, COO uses the 
spherical EPS conditional mass function only in the lower 
mass range [0, Mo/2] to generate the progenitor abundance. 
By construction, then, the shape of the progenitor mass 
function in the upper mass range, [Mo/2, Mo], is symmetric 
with the lower half and fails to match accurately the asym- 
metric EPS 0(M|M o ). 



5 THREE CONSISTENT MONTE CARLO 
ALGORITHMS 

In this section, we present three Monte Carlo algorithms 
that all satisfy the criterion for consistency discussed in 
§3.1| and will therefore accurately reproduce the EPS pro- 
genitor mass function <j)(M\Mo). We introduce the common 
setup for our methods in §5.1| and discuss in detail how each 
method assigns the ensemble of progenitors to descendants 
in g5T2|-[5T4l 

To help the reader follow our discussions, we provide a 
summary of the breakdown of the merger configurations for 
the three new algorithms in Table [3] and the accompanying 
Fig. © 

Although the standard practice in the community has 
been to generate merger trees using the spherical EPS 
model, we emphasize that the Monte Carlo algorithms can 
be applied to the ellipsoidal EPS model as well. In fact, 
since the ellipsoidal model matches the unconditional mass 
function in simulations better than the spherical model, we 
would expect the ellipsoidal EPS to also match better the 
progenitor statistics in simulations. We will therefore present 
our results for both the spherical and ellipsoidal EPS models 
in parallel below. 



5.1 The Common Setup 

5.1.1 Basic Features 

Our Monte Carlo algorithms for growing consistent merger 
trees all share the following implementation framework. We 
begin at redshift and build the merger tree backwards in 



cosmic time. We typically choose a large descendant halo 
mass range (M = [1O 6 M , 1O 15 M ]) and a small simula- 
tion time-step (Az m 0.02 at low 2; see discussion below) to 
achieve a high resolution tree and a large dynamic range in 
the progenitor mass. For a given descendant halo, we first 
compute which mass bin it belongs to, and then obtain its 
progenitors across a single time-step using the distribution of 
merger configurations specific to each algorithm (described 
in the next three subsections). The progenitors then become 
descendants in the next time-step, and this process is re- 
peated to build up the higher tree branches. 

To be specific, a merger configuration here is defined as 
a set of progenitor masses that form a descendant halo of a 
given mass in one time-step. For example, for a descendant 
halo of mass Mo, one merger configurations may include 
only two progenitors of mass O.6M0 and 0.4Mo, while an- 
other may contain three progenitors of mass 0.4Mo, 0.3Mo, 
and 0.2Mo. Note that the sum of the progenitor mass in each 
configuration need not equal the descendant mass, and the 
deficit, AM, is implicitly attributed to sub-resolution pro- 



genitors (see £ 3.3 1 . Different Monte Carlo algorithms have 
different distributions of merger configurations and progen- 
itor multiplicities for each descendant bin. For convenience, 
we call the most massive progenitor in a merger configura- 
tion the primary progenitor, and the rest of the progenitors 
the secondary progenitors. 

Our basic implementation is applicable to both the 
spherical and ellipsoidal EPS models. We find a particu- 
larly efficient choice of time-step to be the one correspond- 
ing to a constant difference in the barrier height Aui(z) — 
uj(z + Az) — lo(z), as is used in ^4]for the four old algorithms. 
For the spherical case, the progenitor mass function eq. |3| 
depends on time only through Aui(z) and is therefore iden- 
tical for all redshifts when the same Auu(z) is used. Thus 
we only have to generate the merger configurations in the 
spherical case across a single time-step once. For the ellip- 
soidal case, however, the progenitor mass function eq. |5| 
not only is a function of Au(z) but also depends explic- 
itly on z. For each Monte Carlo algorithm, it is therefore 
necessary to generate and store the merger configurations 
and their probabilities for both descendant halos of differ- 
ent masses and several redshift bins. In practice, since the 
redshift dependence of eq. |5| is weak, typically fewer than 
~ 20 redshift bins are required. 



5.1.2 Important Progenitor Mass Scales 

A number of natural mass boundaries play critical roles in 
the construction of our algorithms. These mass scales de- 
marcate the regions with different progenitor multiplicities, 
as illustrated in Fig. [6] and discussed in detail in the next 
three subsections. 

(i) The resolution scale M rcs and its complement Mo — 
M res are two obvious boundaries, as is the half descendant 
mass Mo/ 2 discussed in the context of binary mergers in 
Sec. 3.2 We generally choose a small M ros (typically M rea = 



O.OOlA/o) for numerical precision and keep track of all the 
progenitors down to this limit at each time-step, 
(ii) The mass qMq given by 



Mi 



<t>{M\Mo)dM = 1 



(14) 



a M 
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a, Spherical 



ft, Spherical 



, Ellipsoidal 



fi, Ellipsoidal 
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Figure 5. a and as functions of the look-back time Az at red- 
shift zero. The red solid, blue dotted, and black dashed curves 
are for descendant halos of 1O 12 M , 1O 13 M , and 1O 14 M re- 
spectively. The label in each plot indicates the quantity {a or fi) 
shown and the EPS model (spherical or ellipsoidal) used. 



defines the range of progenitor mass over which every de- 
scendant halo is guaranteed to have one progenitor with 
M € [aM ,M ]. Table [2] lists the values of a for both the 
spherical and ellipsoidal progenitor mass functions for three 
descendant masses; a is seen to range from 0.361 to 0.448. 

(iii) The mass fiMo demarcates where the asymmet- 
ric progenitor mass function self-intersects: ^>(/xMo|Mo) = 
4>{Mo — (jlMo\Mo) with fi > 0.5. For binary merger configu- 
rations of the form Mo = M1+M2, 4>(Mi\M ) > 0(M 2 |M o ) 
when Mi < /j,M and </>(M 2 |M ) > 0(Mi|M o ) when Mi > 
/xMrj. This mass scale is illustrated in the pop-up in Fig. [I] 
Table [2] shows that fi w 0.956 to 0.977. 

Fig.[5]shows a and fi as functions of the look-back time 
Az for three descendant halo masses (1O 12 M , 10 13 M o , 
1O 14 M ) at redshift zero. According to the figure, a and 
/x have well defined constant values when Az is less than 
about 0.05, a natural upper limit of time step-size for a 
Monte Carlo simulation to achieve convergence in both the 
spherical and ellipsoidal EPS models. 



5.1.3 Mass Bins 

To help the reader reproduce our Monte Carlo algorithms, 
we discuss our distribution of mass bins. 

We divide the descendant mass range 10 6 ^ M ^ 
1O 15 M into ~ 100 logarithmic descendant bins. Halos that 
fall into the same descendant bin are assumed to have the 
same distribution of single-time-step merger configurations 
that are computed using the central (in logarithmic scale) 
value of the bin as the descendant mass. The progenitor 
masses in a merger configuration are recorded in the form of 
ratios to the descendant halo mass, instead of their absolute 
masses. This allows us to correct for the (small) difference 




M„-M re « M 



Figure 6. A schematic summary of how the three new algorithms 
proposed in this paper assign progenitors to descendants in a sin- 
gle time-step (see §5). The regions are shaded according to the 
progenitor multiplicity (marked by N p — > 1) and the mass ranges. 
See Table [3] for a description of each shaded region and the frac- 
tion of descendants that belongs to each region. The numbers 
quoted in this plot are from the ellipsoidal EPS model. The axes 
are in arbitrary units, though the horizontal axis is drawn to be 
symmetric about Mo / 2 and the vertical axis is assumed to be log- 
arithmic. Important characteristic progenitor masses are labeled 
on the horizontal axis (see §5.1.2 for discussion). The dashed line 
in panel A plots </> sym , the reflection of the left side of </>(M\Mo). 



between the descendant halo in question and the central 
mass of its bin. 

For a given descendant mass Mo, its progenitor mass 
range [M ros ,Mo] is divided into a certain number of mass 
bins to facilitate the process of forming merger configura- 
tions. Interestingly, we note that simply dividing the whole 
progenitor mass range into evenly spaced logarithmic bins 
is not accurate, as discussed in §4.2| This is because the 
simplest logarithmic binning assigns very few bins to the 
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Spherical EPS Ellipsoidal EPS (z=0) 



M (M Q ) 


10 12 


10 13 


10 14 


10 12 


10 13 


10 14 


a 


0.421 


0.448 


0.435 


0.361 


0.384 


0.372 


/" 


0.977 


0.977 


0.970 


0.974 


0.970 


0.956 



Table 2. Values of the progenitor mass scales a and fi discussed 
in §5.1.2| for the spherical and ellipsoidal EPS models for three 
descendant masses (10 12 , 10 13 , and 1O 14 M ) and Az = 0.02, 
where aMo is defined such that tj>(M\Mo) dM = 1 and fiMo 

is defined such that (p(fiM \M ) = <j>(M -fj,M \M ) with n ^ 0.5. 



mass range of [Mo/2, Mo], which requires many mass bins 
to sample accurately the shape of the sharply peaked (at 
around Mo) progenitor mass function for a small time-step. 
To give the peaked region more fine structures, we find a 
simple way: the mass range of [M rcs , Mo/2] is divided into 
evenly spaced logarithmic bins, and its reflection about the 
mid point Mo/ 2 determines the binning on the right side of 
the mid point. Mathematically, it can be stated as follows: 
The progenitor mass range [M rcs , Mo] is divided into 27V + 1 
logarithmic mass bins. The i th (i = 0, 1, 2, 2N) bin spans 
[M i+1 ,M% where M i is defined as follows: 

r M if i = 0; 

M i = < exp[lnM rcs + Ax (2/V + l-i)] ifi>/V + l; 
i M - M 2JV+2 -* if 1 < i < N. 

and A = (ln(M /2) - In M ICS )/N. 

The average number of progenitors (per descendant 
halo) in the i th bin is called TV 1 , which is equal to 

J^i+i <(>(M\M )dM. Note that N' is not an integer. For 
i ^ 1, we choose the mean mass M' of the i th bin to 
be \[WW+ l . The progenitor mass function often changes 
rapidly across the th bin so we do not assign it a mean mass. 
Instead, whenever a progenitor of the th bin is needed, we 
generate a probabilistic progenitor mass according to the 
progenitor mass function inside this bin. 

5.2 Method A 

We first attempt to resolve the asymmetry problem in the 
EPS progenitor mass function 4>{M\Mo) by assuming that 
the primary progenitors in the symmetric part </> S ym in 
eq. (JsJ) are paired up with secondary progenitors to form 
binary mergers such that Mo = Mi + M2. This is done 
so long as the smaller progenitor is above the mass reso- 
lution of the Monte Carlo simulation, i.e. M2 ^ M res and 
M\ Mo — M res . If M2 < M rea , then the second progenitor 
is discarded and Mi is assumed to be a single progenitor 
(the darkest grey region marked 1 — > 1 in Fig. [6] A). The re- 
maining primary progenitors in the asymmetric part aS ym 
are assumed not to pair up, i.e. each descendant halo has a 
single progenitor (the lightest grey region marked 1 — *• 1 in 
Fig. ©A). 

In practice, we generate the merger configurations of a 
descendant halo of mass Mo at each time step by repeating 
these two simple steps: 

(i) Draw the primary progenitor mass Mi from the mass 
range [Mo/2, Mo] of the progenitor mass function. 

(ii) If Mi > Mo — M tes , no more progenitors are gener- 



ated; if Mi ^ Mo — M res , the probability of having a second 
progenitor of mass M2 = Mo — Mi is set to 

= t?Wm(Mi|M ) = 4>(Mo -Mi[Mq) 
' ~ 4>{Mi\M ) ~ 4>{Mi\Mo) ' ( ' 

Then, drawing a random number x between and 1 allows 
us to determine whether a secondary progenitor should be 
generated. If x < r, M2 is assigned as a secondary progeni- 
tor; otherwise Mi is left as the sole progenitor. 

We point out two subtleties with this algorithm. First, 
r is not always ^ 1. It is true that r is below 1 for most 
of the relevant mass range Mi £ [Mq/2, fiMg] (see Fig. [6] A 
and Table [5} since the left side of <^>(M|Mo) is slightly lower 
than the right side. But when Mi > fiMo, we find that 
r > 1, implying that on average more than one secondary 
progenitors should be generated to couple with the primary 
progenitor Mi , and we must generate merger configurations 
with multiple progenitors. To accommodate this feature, for 
each Mi that satisfies Mi G [fiMo,Mo — M res ], we gener- 
at^jeither int(r) or int(r) + 1 secondary progenitors of mass 
Mo — Mi according to whether a random number between 
and 1 is larger or smaller than r — int(r). Note that the re- 
sulting merger configurations do not conserve mass exactly 
because the sum of the progenitor masses is slightly larger 
than the descendant mass. Typically most of these config- 
urations only end up with 3 or 4 progenitors as r <J 2 for 
Mi < 0.999M and Az < 0.02. 

The second subtlety with method A is that since the 
total number of progenitors in the mass range of [Mo/2, Mo] 
(which is equal to J^°^ 2 (j)(M\Mo)dM) is always slightly 
smaller than one (typically by 0.2% to 0.4% for Az = 0.02; 
recall from Table[5]that aMo < Mo/2), it is possible that we 
sometimes cannot assign any progenitors to a given descen- 
dant halo. When this happens, the descendant halo does not 
have any progenitor halos and is a "0 — > 1 event". 

For a thorough description of our algorithm A, we list 
below all the possible merger configurations and their fre- 
quencies of occurrence for descendant halos at z = over 
a single simulation time-step Az — 0.02 and mass resolu- 
tion M ros = O.OOlAfo- This information is also summarized 
in Table [3] and Fig. [6] In general, the relative frequencies of 
different merger configurations are insensitive to the descen- 
dant mass Mo but do depend on the Az and M ros used in the 
Monte Carlo simulation. For example, the fraction of 1 — > 1 
events increases as Az decreases; and if M res is chosen to be 
larger than (1 — /i)Mo ~ 0.03Mo, then there are no 3 — > 1 
or 4 — > 1 mergers at each time-step and mass conservation 
is exactly respected. 

I. About 12% in the ellipsoidal model (21% for spher- 
ical) of descendant halos have two progenitors each. These 
are binary pairs drawn from the symmetric part of the pro- 
genitor mass function ci S y m , where Mi 6 [Mo/2, pMo] and 
M = Mi + M 2 (Fig. [6] A ■). 

II. About 69% (60%) of descendant halos have only one 
progenitor each. The majority ( J> 99%) of these descendants 
originally have binary progenitors but the smaller progen- 
itor is discarded since M2 < M rcs (i.e. Mi > Mo — M res ) 
(Fig. [6] A ■). The rest ( <, 1%) of these descendant have 

5 Here int(r) is defined to be the largest integer n that satisfies 
n ^ r. 
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Figure 7. Comparison of the progenitor (or conditional) mass functions <fi(M, z\Mo, 20) generated using the three new Monte Carlo 
algorithms introduced in |J5] A (red solid), B (green dashed), and C (blue dotted), and the predictions of the spherical EPS model (black 
solid). The four panels show four look-back redshifts (z — zo = 0.24,2.07,7,15) for a descendant halo of mass 1O 13 M at z = 0. For 
clarity, we plot in the sub-panel below each panel the ratios of the Monte Carlo result and the EPS prediction. All three algorithms are 
seen to match very closely the spherical EPS prediction at all redshifts. At z = 7 and 15, the slight underestimates of the progenitor 
abundances at M/Mq < 10~ 4 are primarily due to the fact that we trace a halo's progenitors only down to 0.001 of the halo mass in 
each small time-step in our Monte Carlo simulations. 




0.01 0.1 1 10- 5 10- 
M/M (M = 10 13 M o ) 

Figure 8. Same as Fig.[7]except both the Monte Carlo and analytic results are now generated from the ellipsoidal instead of the standard 
spherical EPS model. The Monte Carlo methods use eq. (JsJ as the progenitor mass function for each time step. The analytic results 
are calculated using the integral equation proposed by Zhang & Hui (2006). The agreement is again excellent, indicating that our new 
Monte Carlo algorithms work well in reproducing the EPS progenitor mass function regardless if the EPS model is based on constant (i.e. 
spherical collapse) or moving barrier (ellipsoidal) random walks. For completeness, we include the results from the ellipsoidal version of 
the KW93 method (orange dash-dotted). At z=0.24, the slight progenitor overabundance at the low mass end is due to the approximate 
nature of eq.[5] At z=7 and 15, the slight underestimates of the progenitor abundances are due to both the mass resolution issue as stated 
in the caption of Fig. [7] and the barrier intersection problem of the ellipsoidal EPS model, which prevents us from tracing progenitors 
that are much smaller than the typical halo mass of the same redshift. 
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Method 


N p 




% Desc. 
(spher.) 


% Desc. 
(ellip.) 


Description 


Kev 


A 


-+ 


1 


0.3% 


0.4% 


Descendants with no progenitors because / <j> dM < 1 

J jV/q / 2 


N/A 




1 -> 


1 


60% 


69% 


Mi G [Mo — M res ,Mo]: binary-turncd-singles due to M2 < M res 


A ■ 




1 -> 


1 


0.8% 


0.4% 


Mi S [Mo/2,/xMo]: S ym < results in unpaired primary progenitors: 

AM > M res 


A □ 




2 


1 


21% 


12% 


Mi £ [Mo/2, fiMo] and Mo = Mi +M2: binary pairs generated from </> S ym 


AH 




3+- 


• 1 


18% 


18% 


Mi £ [aiMo, Mo — M res ]: sym > results in excess secondary progenitors. 
M < Mi + M 2 +M 3 + ... 


A ■ 


B 


1 -> 


1 


60% 


69% 


Mi 6 [Mo — Mres j Mo]: binary-turned-singles due to M 2 < M rea 


B ■ 




2 -» 


1 


20% 


14% 


Binary paired progenitors generated by the iterative algorithm of §5.3| 
M ^ Mi + M 2 


B ■ 




3+ - 


» 1 


20% 


17% 


Mi G [/xMo,Mo — M res ]: identical to 3H — > 1 configuration in method A 


B ■ 


C 


1 -> 


1 


60% 


69% 


Mi e [M - M res , Mo]: binary-turned-singles due to A/ 2 < M TCB 


C ■ 




1 -> 


1 


35% 


29% 


All secondary progenitors have already been assigned to smaller primary 
progenitors: these remaining primary progenitors have AM > M re s 


c ■ 




3 -* 


1 


0.1% 


0.01% 


Merger configurations with 3 progenitors 


CM 




4 -> 


1 


2% 


0.3% 


Merger configurations with 4 progenitors 


c n 




5+ - 


> 1 


2.9% 


1.7% 


Merger configurations with 5 or more progenitors 


c □ 


KW93 


1 -» 


1 


60% 


69% 


Mi e [M - Mres, Mo]: binary-turned-singles due to Af 2 < Afros 


N/A 




1 -» 


1 


15% 


9% 


Merger configurations with a single progenitor with Mi < Mo — M res 


N/A 




2 -► 


1 


11% 


9% 


Merger configurations with 2 progenitors 


N/A 




3+ - 


> 1 


14% 


13% 


Merger configurations with 3 or more progenitors 


N/A 



Table 3. A summary of our three new Monte Carlo methods discussed in ij5]and the method of KW93. The percentages indicate the 
fractions of descendants with N p progenitors in a given method, computed for Mo = 10 13 Mq and M re s = O.OOlMo for a single time-step 
Az = 0.02 in both the spherical and ellipsoidal EPS models. They are representative of the merger configuration distributions for other 
descendant halo masses Mq. 



progenitors with Mi £ [Mo/2, fiMo] and originate from the 
small asymmetric part </> a sym of the progenitor mass function 
where r < 1 (Fig. [6] A □). 

III. About 18% of descendant halos have three or four 
progenitors each, typically consisting of one massive pro- 
genitor and two or three very small secondary progenitors 
( < (1 - n)M ~ 0.03Mo). These have Mi £ [fiM ,M - 
Afres] (Fig. [6] AH). 

IV. About 0.4% (0.3%) of the descendants have no pro- 
genitors due to the sharp cutoff of the primary progenitor 
mass at Mo/2 discussed above. 

The red solid curves in Fig. [7] compare the progeni- 
tor mass functions from this Monte Carlo algorithm with 
the analytic eq. (J3J) of the spherical EPS model. Fig. [8] 
shows the same thing except everything is for the ellip- 
soidal EPS model, where we use eq. |5| to compute the pro- 
genitor mass function for each small simulation time-step. 
Both figures show excellent agreement (< 10% deviation) at 
z\ — zo = 0.24, 2.07, 7, and 15 for a descendant halo of mass 
1O 13 M0 at zo — 0. We have tested other descendant masses 
(1O 12 M < M < 1O 14 M ) and found equally good agree- 
ment. This agreement also provides numerical verification of 
the criterion introduced in §3.1| 



5.3 Method B 

Two features in method A may seem unnatural. First, as 
shown in Table [3] and discussed in the previous section, a 
small fraction (~ 0.3% to 0.4%) of the descendant halos 
in method A are not assigned any progenitors in one time- 
step because J™° 4>(M\M )dM w 0.997 (for Az = 0.02 
and a large range of Mo) and is not exactly unity. It is 
important to note that though these descendants are rare, 
one cannot remove them from method A by modifying the 
normalization of <f>(M\Mo) in the mass range of [Mo/2, Mo], 
as such a modification is amplified with iterations and leads 
to a large error in <j>(M\Mo) after many time-steps. 

Second, due to the asymmetry in the EPS <f>(M\Mo), we 
have assigned a small fraction (0.4% to 0.8% for parameters 
used in Table [3| of the descendant halos to 1 — ► 1 events. 
There is therefore a small chance that a progenitor of mass 
comparable to half of the descendant mass does not have 
any companions, corresponding to a large deficit between 
the mass of the descendant halo and the total mass of its 
progenitors. 

The first feature can be avoided by decreasing the lower 
limit of the mass range from which the primary progenitor is 
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drawn from Mq /2 to aMo , where a is defined in eq. ( f 4 1 and 



ranges from a ~ 0.36 to 0.45 in Table [2] The second feature 
can be altered by distributing the secondary progenitors in 
a slightly different way. These options motivate us to invent 
Method B with the following set up: 

1. We assume the primary progenitor mass lies in the 
mass range [aMo, Mo]. This condition guarantees that every 
descendant halo has a primary progenitor of mass > aMo 
due to the definition of a. 

2. We then assign secondary progenitors to primary 
progenitors from the left side of aMo- For simplicity, when- 
ever possible, we make only binary configurations, each of 
which contains one primary and one secondary progenitor. 
We start with the primary and secondary progenitor bins 
that share the aMo boundary (i.e. nearly equal-mass pairs) 
and work our way outwards to the Mi 3> M2 pairs. This 
is a natural decision as this way of pairing the primary and 
secondary masses minimizes the difference between Mo and 
Mi + M2 . Specifically, for a given Mi bin, we determine its 
binary companion's mass M2 from 

<f>(M\M ) dM = / cj>(M\M )dM, (16) 

' M-2 J aM 

which guarantees that we always have an equal number of 
secondary progenitors to pair up with the primary halos. 
Note that since a < 0.5 it is generally true that Mo > 
Mi + M 2 . 

3. One caveat with step 2 above is that this simple 
binary pairing scheme works for a large range of masses but 
needs to be modified near the end points when Mi is close 
to Mo and M2 <C Mi. This is because the scheme starts 
out with nearly equal-mass pairs at Mi ~ M2 ~ aMo and 
Mi +M2 < Mo, and the asymmetric shape of the progenitor 
mass function is such that the method produces pairs with 
increasing Mi + M2 as we move outward from aMo- The 
equality Mi + M2 = Mo is reached when Mi is slightly 
larger than /iM (typically at 0.99Mq), beyond which there 
are more secondary progenitors left to be paired than the 
primary ones. We therefore stop the binary pairing when 
AIi + M2 = Mo is reached. From this point on, we instead 
use the same multiple merger configurations as in method 
A. For simplicity in the following few paragraphs, we denote 
this transitional Mi as //Mo. 

In summary, methods A and B are closely related and 
are compared side- by-side in Table [3] and Fig. [6] They have 
identical merger configurations in the following regions: 

I. The high- Mi region Mi £ [Mo - M rcs ,M ], where 
60% to 70% of descendant halos belong. The secondary pro- 
genitor is below M rcs , so Mi is effectively the sole progenitor 
(i.e. N p — 1) for these descendants 

II. The region Mi G [ii 1 M , M ~M ICS ] (fj,' replaced by fi 
in method A), where 17% to 20% of descendant halos belong. 
These descendants all each have 3 or more progenitors (N p = 
3+). 

Methods A and B differ in the following regions: 

III. The binary pairing algorithm used in method B 
removes the sliver of 1 — ► 1 configurations in the Mi £ 
[Mo/2, fj.Mo] region in method A (A □) and redistributes the 
binary merger configurations in this region (A fl) to yield a 
robust set of binary configurations between aMo ^ Mi ^ 
fx Mo (B □). This affects ~ 20% of the descendant halos. 

IV. Since the primary progenitor mass range extends 



down to aMo instead of Mo/2, method B does not have any 
of the — > 1 configurations that are present in method A. 

The green dashed curves in Figs. [7] and [8] compare the 
progenitor mass functions from method B with the analytic 
predictions of the spherical and ellipsoidal EPS models, re- 
spectively. The agreement is again excellent (< 10% devia- 
tion) at zi — zo — 0.24, 2.07, 7, and 15 for a descendant halo 
of mass 1O 13 M at z = 0. 

Finally, we note that mass is not strictly conserved 
for the multiple merger configurations generated in the 
Mi e [liM ,M - Mres] region of method A and the Mi G 
[n'M ,M - M tea ] region of method B (Fig. [6] A ■, B ■). 
These configurations have more than one companion of mass 
Mo — Mi , making the total mass of the progenitors slightly 
above the descendant halo mass. This issue is due to the 
rapid rise of the progenitor number as the secondary progen- 
itor mass approaches zero. In principle, the small progenitors 
( < (1 — ;u) Mo) that are causing this problem can be re- 
distributed and combined, e.g. , with progenitors in some of 
the 1 — * 1 and 0—^1 merger configurations in method A, or 
with some binary configurations of total masses smaller than 
the descendant mass in method B, to form multiple merger 
configurations that obey mass conservation (this, in fact, 
is what happens in method C below, where mass conserva- 
tion is strictly respected). We have checked that this can be 
done successfully without violating mass conservation down 
to very small M res and find that in practice, these modifica- 
tions do not introduce significant changes to the statistical 
properties of the halo merger histories. We have therefore 
chosen to present the simpler version of each model. It is 
also worth noting that in the EPS theory, mass conserva- 
tion only has to be obeyed statistically and is not required 
for individual merger configurations. 



5.4 Method C (Multiple Mergers) 

As shown in Table [3] methods A and B both produce com- 
parable number of descendants with binary (N p = 2) and 
multiple (N p = 3+) progenitors in a single time-step. The 
importance of multiple merger configurations have been em- 



phasized by a number of authors (e.g., Kauffmann & White 
[Somerville fc Kolatt||1999| [Neistein fc Dekel||2008bT 



1993 



It is therefore interesting to explore the relative importance 
of binary vs multiple mergers by relaxing the binary as- 
sumption. Our method C is designed for this purpose. More 
specifically, this method does not have any restrictions on 
the number of progenitors in each merger configuration. We 
only require that the total progenitor mass of every merger 
configuration be smaller than (or equal to) the descendant 
halo mass. 

We now describe method C: 

1. To prevent the formation of — > 1 merger configura- 
tions we mimic the setup of method B and choose to draw 
primary progenitors from the mass range Mi G [aMo, Mo]. 
Thus methods B and C share the same distribution of pri- 
mary and secondary progenitor mass bins. 

2. As with method B, we form merger configurations by 
assigning secondary progenitors to progenitors in primary 
bins. Every primary bin starts with one merger configura- 
tion: that which contains only the primary progenitor itself, 
and has a probability iV con f equal to the number of primary 
progenitors in the bin. The assignment of secondary progeni- 
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tors to primary bins is done in order of decreasing secondary 
progenitor mass. For each secondary bin, we scan the pri- 
mary bins in order of increasing primary progenitor mass to 
find configurations with room to hold at least one secondary 
progenitor from the bin in question (recall that we require 
the sum of progenitor masses to never exceed the descendant 
mass). 

3. When a valid configuration is found, we always as- 
sign the maximal number of secondary progenitors to that 
configuration. For example, suppose we start to assign sec- 
ondary progenitors from a bin with central mass M2 (say 
there are N2 such progenitors in this bin), and find a valid 
configuration of probability N conz and total progenitor mass 
Mtot- The maximum number n ma x of secondary progenitors 
that can be added into each realization of this configuration 
is equal to int[(Mo — Mtot)/M<i\. Therefore, we can maxi- 
mally assign N max = n max x N CO nf secondary progenitors to 
this configuration. 

I. If N m ax > N2, we break the configuration into two: 
one contains the original set of progenitors, with a probabil- 
ity equal to (1 — Ni/Nmax) x N con r, the other contains the 
original set of progenitors plus n max secondary progenitors 
of mass M2, with a probability equal to (N2/N max ) x N conl . 
In this case all the secondary progenitors of the current sec- 
ondary bin are assigned. 

II. If N max ^ N2 we simply add the Umax secondary 
progenitors of mass M2 to the configuration, and update the 
list of progenitors in the configuration. N con { , the configura- 
tion's probability does not change. The number of remaining 
secondary progenitors to be matched is now N2 — N m ax, and 
we continue our search across merger configurations (in or- 
der of increasing primary progenitor mass) until all of them 
have been assigned. 

Once a secondary bin is fully assigned, we move on to 
the next secondary bin (of a slightly smaller mass) and re- 
peat the same assignment procedure. As this process goes on 
all configurations are gradually filled with secondary progen- 
itors of smaller and smaller mass. For technical convenience, 
the number of configurations in each primary bin and the 
number of unique progenitor masses in each configuration 
are both limited to be fewer than 6. In practice, we find that 
this setup allows us to successfully assign all secondary pro- 
genitors in the mass range [M rcs , aMo], even when the mass 
resolution of each time step is as low as M rcs = O.OOlMo. 

In fact this dense packing of secondary progenitors into 
primary bin configurations manages to distribute efficiently 
all secondary progenitors in [M res , aMo] in only a fraction of 
the available primary progenitors. As seen in Fig. [6]C, only 
2% (5% for spherical) of the primary progenitors (at the low 
mass end) are grouped with secondary progenitors and the 
remaining 98% (95%) are 1 — > 1 events. We note that even 
though there are far more secondary progenitors than pri- 
mary progenitors, this is possible because many secondary 
progenitors have exceedingly small masses and can be effi- 
ciently distributed into the mass reservoirs of relatively few 
primary progenitors. 

The execution of method C is as follows: 

(i) Generate a primary progenitor Mi from the mass 
range [aMo, Mo] of the EPS progenitor mass function. De- 
termine which primary bin contains Mi. 

(ii) If Mi > Mo — Mrcs, no more progenitors are gen- 
erated; if Mi ^ Mo — Af rcs , a random number determines 



which merger configuration to choose according to the prob- 
ability distribution of all possible configurations associated 
with the given primary bin. The progenitors of the chosen 
configuration are then generated. 

For a better understanding of method C, we show in 
Table [3] and discuss below all the possible merger configu- 
rations and their frequencies of occurrence for descendant 
halos (regardless of their masses) at z — 0, assuming time- 
step Az — 0.02 and mass resolution M rcs = O.OOlAfo: 

I. About 98% (95% for spherical) of the descendant ha- 
los have only one progenitor each. 

A) About 2/3 of these descendants' progenitors are 
within the resolution limit of the descendant mass (i.e. 
Mi > M - M-os, see figure [6] C ■). 

B) The remaining 1/3 of these descendant halos' pro- 
genitors have masses below Mo — M rcs . As discussed above, 
these massive primary progenitors are not assigned any sec- 
ondary companions because all the available secondary pro- 
genitors are maximally assigned to the less massive primary 
bins. Note that this region extends to masses below [iM^ ( 

cm). 

II. For the remaining primary progenitor bins, there are 
no configurations having only two progenitors. All in all, 
0.01% (0.1% for spherical) of all descendants have three pro- 
genitors (C ■); 0.3% (2%) have four progenitors (C ■); 1.7% 
(2.9%) have five or more progenitors (C □). The progenitor 
count for a given configuration can be rather large reaching 
values of more than 100. 

As in methods A and B, the values quoted above de- 
pend on Az and M rcs . They also depend on the maximal 
number of configurations allowed in each primary bin and 
the maximal number of unique progenitor masses allowed in 
each configuration. 

The blue dotted curves in Figs. [7] and [5] compare the 
progenitor mass functions from this Monte Carlo algorithm 
with the analytic predictions of the spherical and ellipsoidal 
EPS models, respectively. They again show excellent agree- 
ment (< 10% deviation) at 21 - z = 0.24, 2.07, 7, and 15 
for a descendant halo of mass 10 13 M Q at z = 0. 



6 COMPARISON OF HIGHER-MOMENT 
STATISTICS IN ALGORITHMS A, B, C, 
AND KW93 

We have designed Monte Carlo algorithms A, B, and C for 
constructing merger trees that can accurately reproduce the 
EPS prediction for the progenitor mass function <p(M\Mo) 
across each individual time-step. According to the discus- 
sion in §3.1[ these methods should then accurately generate 
the progenitor mass function at any look-back time in any 
number of time-steps. Figs. [7] and [5] show that this is indeed 
the case for both the spherical and ellipsoidal EPS models. 
Including KW93, there are now four methods that are com- 
pletely consistent with the EPS 0(M|M o ). The results of 
the ellipsoidal version of KW93 have been shown in Fig. [8] 
as well. 

Despite this agreement, we recall that the progeni- 
tor mass function is only one of many statistical prop- 
erties of a halo merger tree. Even though all four algo- 
rithms are degenerate in <f>(M\Mo), they are likely to (and 
should) differ in their predictions for other statistical quan- 
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tities. Here we investigate two such quantities as an illus- 
tration: (i) cp fNp ^ (M\Mo), the progenitor mass function for 
the subset of descendant halos that have N p progenitors. 
The sum of (Jv f ) (M\M ) over all N p is equal to 0(M|M o ). 
(ii) <f)^ Hh) (M\Mq), the distribution of the ith most mas- 
sive progenitor of each descendant halo. Again, the sum of 
4> {Hh) (M\M ) over all i is equal to cj>(M\M ). These two 
statistics are two obvious ways of decomposing the total 
<j)(M\Mo) into individual moments: (p^ Np ' separates flour- 
ishing trees from quiescent trees, while <j)^ th ^ compares the 
individual distributions of the primary, secondary and more 
minor progenitors, which are relevant for modeling galaxy 



formation through mergers (see also Parkinson, Cole & Helly 
2008). Other statistics such as the distributions of halo for- 



mation time and last major merger time (e.g 
Cole & Helly 



2008 



Parkinson, 

Cole et al.|2008"l|Moreno fc Sheth|20"07] ) 



and the factorial moments of the partition function (e.g., 
Sheth & Pitman 1997 Sheth & Lemson||1999 1 are also use- 



ful. Some of these will be examined in our next paper. 

To compute these moments, we set the descendant halo 
at redshift zero to be 1O 13 M0, and the mass resolution to be 
4 x 1O 1O M0. The results are plotted at two look-back times 
(zi-zo = 0.51, 2.07) in Figs. [9p2| where Figs.[9]and[To]show 
(JVp) (M|M o ) for the spherical and ellipsoidal EPS models, 



respectively, while Figs. 



11 



and 



12 



show (i * fe) (M|Mo). In 



each figure, results from our three methods (red solid for 
A, green dashed for B, blue dotted for C) and from our 
implementation of KW93 (orange dash-dotted) are shown 
for comparison. These figures clearly indicate that methods 
A, B, C, and KW93 generate distinct predictions for these 
specific moments of the progenitor mass distribution. Some 
of the notable differences are: 

1. Method C produces a much lower amplitude for the 
N p — 2 and 3 moments than methods A and B. This is 
because C is designed to be a multiple-merger algorithm that 
effectively does not generate any binary configuration in one 
individual time-step (note the absence of the N p = 2 entry 
for method C in Table |3j. This feature can been seen by 
the absence of blue short-dashed curves in the N p = 2 and 3 
panels in Figs. [9| and |10| i.e., there are almost no descendant 
halos having only two or three progenitors in method C at 
z = 0.51. By contrast, methods A and B have a wealth of 
descendants with binary progenitors at these redshifts. 

2. The removal of the binary assumption in method C 
leads to many features in the moments of the progenitor 
distributions. By contrast, the predictions from A and B 
are mostly power-laws, or at least smooth functions, in the 
progenitor mass. This difference is due to the fact that the 
merger configurations in the binary methods are much more 
regulated than those in the non-binary method: a binary 
configuration contains only two progenitors, the total mass 
of which is always quite close (if not equal) to the descen- 
dant mass, whereas the distribution of progenitor masses in 
a multiple configuration can have various forms, which can 
easily affect, e.g. , the ranking of the progenitor masses and 
the number of progenitors. It is interesting to note that the 
predictions of KW93 are fairly smooth functions in spite of 
the fact that it does not assume binary. This is likely be- 
cause the way progenitors are assigned in KW93 effectively 
suppresses the probability of mergers involving multiple pro- 
genitors. 

3. The differences between method A and B are more 



subtle because they are both mostly binary methods. The 
main feature that distinguishes A from B is in the distribu- 
tion of the most massive progenitor (i.e. ith ~ 1) shown in 
the first columns of Figs. |ll| and |12| At the high mass end, 
method B has a slightly broader shape for the primary pro- 
genitor mass than method A. This is expected, because it is 
the case across every time step by construction (the primary 
bins in method B extend down to aMo as opposed to Mo/2 
for method A). At the low mass end, however, there is a 
long tail in the distribution of primary progenitor masses in 
method A, which is not present in other methods. This tail is 
caused by the fact that in method A, there is a small chance 
(~ 0.3%) at every time-step that a primary progenitor com- 
pletely disappears, transferring the rank of "primary" to 
one of the much smaller secondary progenitors. Over sev- 
eral time-steps this rare occurrence affects more and more 
branches of the merger tree and can significantly modify the 
primary progenitor statistics. 

In summary, we have constructed three Monte Carlo al- 
gorithms that can all reproduce closely the progenitor mass 
function of the EPS model (both spherical and ellipsoidal). 
The methods, however, produce significantly different higher 
moments of the progenitor distributions. They are also very 
different from KW93. Either a theoretical model more com- 
plete than the EPS or direct iV-body results will be needed 
to determine which, if any, of the thus-far successful algo- 
rithms is the winner. We will turn to this subject in the next 
paper (Zhang, Fakhouri & Ma, in preparation). 



7 CONCLUSIONS AND DISCUSSION 

Monte Carlo algorithms based on the spherical EPS model 
have been an essential tool for many studies of galaxy and 
structure formation. These algorithms allow one to generate 
realizations of actual halo merger histories starting from a 
limited set of statistical information about dark matter halo 
properties provided by the EPS model. Since the EPS model 
does not uniquely determine many statistical quantities of 
halo mergers beyond the progenitor mass function, there 
is considerable freedom in how to combine progenitors to 
form descendant halos in each time step in a Monte Carlo 
algorithm. 

The emphasis of this paper is on elucidating and quan- 
tifying the ability of a Monte Carlo algorithm to construct 
merger trees that match the analytic progenitor mass func- 
tion of the EPS model (both the spherical and ellipsoidal 
versions). Four main conclusions can be drawn: 

1. We have shown rigorously that to match the EPS 
progenitor mass function accurately at any look-back time, 
it is necessary and sufficient for a Monte Carlo algorithm to 
reproduce the exact progenitor mass function at each time 
step. 

2. We have reviewed and compared the four most fre- 
quently used Monte Carlo algorithms based on the spheri- 
cal EPS model in the literature: Lacey & Cole 1993, Kauff- 
mann fc White|19931|Somerville &: Kolatt|1999| and|Cole"et 



alJ2000| As seen in Figs. [2jiJ all but KW93 only approxi- 
mately reproduce the spherical EPS progenitor mass func- 
tion at each time step, resulting in large deviations from the 
spherical EPS predictions after the accumulation of small 
errors over many time steps. 
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z: 0.51 -> 
0.01 0.1 



0.01 0.1 1 0.01 0.1 1 

M/M (M = W 13 M e Spherical) 

Figure 9. Predictions of algorithms A (red solid), B (green 
dashed), C (blue dotted), and KW93 (orange dash-dotted) for 
(j>( N p}(M, z\Mq, 20), the mass function of progenitors for descen- 
dant halos that have a total of N p progenitors. Two look-back 
redshifts are shown: z — zq = 0.51 (left) and 2.07 (right). For 
each redshift, four representative values of N p are shown (from 
top down). The simulations are for the spherical EPS model and 
assume a descendant halo mass of 10 13 M© at 2 = and mass 
resolution of M ros = 4 X 10 10 M ra . 



Their problems (see Table [T] for details) can be summa- 
rized as: (i) SK99 generally over-estimates the abundances 
of small progenitors by about a factor of two; (ii) LC93 over- 
produces progenitors by a factor of a few when the look-back 
time is large (Az 3> 1); (iii) COO under-predicts the progen- 
itor abundance at the high mass end when the look-back 
time is large. The origin of these discrepancies frequently 
comes from the incompatibility between the binary merger 
assumption used in the Monte Carlo algorithm (e.g. LC93, 
COO) and the asymmetric progenitor mass function of the 
EPS model. 

3. We have designed three new Monte Carlo algorithms 
that all reproduce closely the EPS progenitor mass function 
over a broad range of redshift (z\ — zq up to at least 15) 
and halo mass. Our methods A and B assign binary pairs 
to the symmetric part of <j>(M\Mo) and non-binaries to the 
asymmetric part; the two differ in the mass ranges for the 
most massive progenitors. Our method C, on the other hand, 
completely relaxes the binary merger assumption. The algo- 
rithms are tested for both the spherical and ellipsoidal EPS 
models and the results are shown in Figsj7] and [8] We see 




0.01 0.1 1 0.01 0.1 

M/M Q (M = 10 13 M S Ellipsoidal) 

Figure 10. Same as Fig. [9] except the Monte Carlo results are 
generated from the ellipsoidal instead of the standard spherical 
EPS model. 



that all three methods perform equally well at reproducing 
the respective progenitor mass function at higher redshifts, 
regardless of whether the spherical progenitor mass function 
eq. (J3j) or ellipsoidal progenitor mass function eq. (JsJ is used 
as input. 

4. As emphasized throughout the paper, the EPS model 
only provides a partial statistical description of dark matter 
halo properties; it does not tell us explicitly how to group 
progenitors into descendants in a Monte Carlo realization. 
Therefore, there are different ways to combine progenitors 
into descendant halos in consistent Monte Carlo algorithms. 

We have used our three new algorithms to illustrate 
this exact point. Despite their success in generating merger 
trees that accurately reproduce the EPS progenitor mass 
function, Figs. |9|12| show that the three algorithms make 
significantly different predictions for quantities such as the 
distribution of the most (or the 2 n£ j or 3 r d most) massive 
progenitor masses, and the mass function of progenitors in 
descendant halos with N p (= 1, 2, 3...) progenitors. A theory 
more complete than EPS would be needed to predict these 
higher-order merger statistics and break the degeneracies 
in the progenitor mass function. Alternatively, comparisons 
with A-body simulations should determine which, if any, 
of the three new algorithms is viable. We view the EPS 
models (spherical or ellipsoidal), Monte Carlo algorithms, 
and A-body simulations as three major components in the 
general study of the formation, growth, and clustering of 



© 2006 RAS, MNRAS 000, [l]-[20l 



Growing Healthy Merger Trees 19 



z: 0.51 -> 
0.01 0.1 1 



z: 2.07 -> 
0.01 0.1 



z: 0.51 -t 
0.01 0.1 1 01 



] 00 
10 
1 

0.1 

0.01 

1(1-' 

100 
10 

I 

1 
0.01 

io- 3 

100 
10 

£ i 

0,1 
U.01 

io- 3 

100 
10 

1 

0.1 
0.01 

io- 3 



-Q- 
X 



- 1 1 1- 

r — Method A L = l 1 
: -- Method B 

Method C 1 

- 4 ~ 
i it J 


- i i i- 

! i»=l | 


h+Httj — i 1 1 mi — i 1 1 iiii- 
r j,i, = 2 n 


H+Httj 1 1 Mill 1 1 Mill! 

r % = 2 i 


H+Httj 1 1 1 Nil 1 1 Mill- 

i s a j 


^hhh] — 1 1 Mill 1 1 llllll 


H+Htt| 1 1 1 INI 1 1 1 llll- 

r" hh - s i 

■ i- 


: 1 1 1 in] — i i inn — i i iiiii: 



1 00 




1 00 


1 




1 








1 




1 


0.01 




0.01 


IO -3 




IO" 3 


100 




100 


10 




1 


] 




"1 


1 




1 


01 




01 


i n-3 

1U 




1U 


100 




"1 00 


1 (] 


X 


1 




> = 




1 




n i 

U. 1 


01 




01 


in -3 

1U 






LOO 




100 


10 




10 


] 




1 


0.1 




"1 


0.01 




0.01 


IO" 3 




IO" 3 



0.01 0.1 1 0.01 0.1 1 

M/M (M = 1O 13 M Spherical) 

Figure 11. Same as Fig. [9]except for a different progenitor statis- 
tic: <f>( Hh ^(M, z\Mo, zo), the mass function of the i t h most massive 
progenitor. 



dark matter halos. In this paper we have focused on the first 
two areas, comparing various Monte Carlo algorithms for 
generating halo merger trees and quantifying their abilities 
to consistently match the analytical EPS progenitor mass 
functions over a broad range of mass and redshift. In our 
next paper (Zhang, Fakhouri, Ma 2008b), we will turn to 
comparisons with the Millennium simulation. 

Several recent papers have investigated other Monte 
Carlo methods (see, e.g ., |Parkinson, Cole fe Helly | 2008 
Neistein fc Dekel||2fj)8a| |Moreno fc Sheth||2007| |Neistein 



Dekcl 2008b Although a complete review of these meth- 
ods is beyond the scope of this paper, it is worth pointing 
out some of their features. The method of lMoreno fc Shethl 



( 2007 1 is essentially equivalent to LC93 but is based on the 
ellipsoidal collapse modeQand is discretized in mass instead 
of redshift. The two progenitor masses for each time step are 
assigned using computer-generated random walks and mov- 
ing barriers. Since the asymmetry problem of the progenitor 
mass function is also present in the ellipsoidal model, this 
method does not accurately reproduce the theory-predicted 
progenitor mass function at each time step. Such a dis- 



They use a square-root approximation for the moving barrier 
form, which avoids the barrier crossing problem at the low mass 
end. 




0.01 0.1 1 0.01 

M/M (M = 1O 13 M Ellipsoidal) 

Figure 12. Same as Fig. |11| except the Monte Carlo results are 
generated from the ellipsoidal instead of the spherical EPS model. 



crepancy is amplified with increasing redshift and is indeed 
sho wn in Fig. 5, 6, and 7 of |Moreno fc Sheth| p007| . 



Neistein & Dekel ( 2008b I have proposed a method 



that exactly reproduces the progenitor mass function of the 
spherical EPS model at each time step. This feature alone 
guarantees it to be consistent with EPS at any look-back 
time according to our discussion in §3.1. However, since the 
method requires solving several differential equations with 
nontrivial boundary conditions for the progenitor masses, it 
is technically harder to implement it. 

Finally, the methods described i n |ParMnson, Cole fc| 
Helly | (|2008[) and |Neistein fc Dekel| (|2008a[) are proposed 



to mimic N-body results. They are based on fitting to N- 
body data rather than the EPS theory. It will be interesting 
to compare the predictions for the various merger statistics 
discussed in this paper from these methods with those from 
our ellipsoidal EPS-based methods and from N-body simu- 
lations. This will be done in the next paper. 
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